home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / clang / c1.zip / TAN.C < prev    next >
Text File  |  1987-06-18  |  2KB  |  89 lines
  1.  
  2. /***********************************************************
  3.  *               The TULSA IBM C BOARD                     *
  4.  *                   918-664-8737                          *
  5.  *             300/1200 XMODEM, 24 Hours                   *
  6.  **********************************************************/
  7.  
  8.  
  9.  
  10. #include "math.h"
  11. #include "errno.h"
  12.  
  13. extern int errno;
  14.  
  15. static double tansub();
  16.  
  17. double cotan(x)
  18. double x;
  19. {
  20.         double y;
  21.  
  22.         y = fabs(x);
  23.         if (y < 1.91e-152) {
  24.                 errno = ERANGE;
  25.                 if (x < 0.0)
  26.                         return -HUGE;
  27.                 else
  28.                         return HUGE;
  29.         }
  30.         return tansub(x,y,2);
  31. }
  32.  
  33. double tan(x)
  34. double x;
  35. {
  36.         return tansub(x, fabs(x), 0);
  37. }
  38.  
  39. #define P1 -0.13338350006421960681e+0
  40. #define P2 +0.34248878235890589960e-2
  41. #define P3 -0.17861707342254426711e-4
  42. #define Q0 +1.0
  43. #define Q1 -0.46671683339755294240e+0
  44. #define Q2 +0.25663832289440112864e-1
  45. #define Q3 -0.31181531907010027307e-3
  46. #define Q4 +0.49819433993786512270e-6
  47.  
  48. #define P(f,g) (((P3*g P2)*g P1)*g*f + f)
  49. #define Q(g) ((((Q4*g Q3)*g Q2)*g Q1)*g Q0)
  50.  
  51. #define YMAX 6.74652e09
  52.  
  53. static double tansub(x, y, flag)
  54. double x,y;
  55. {
  56.         double f, g, xn;
  57.         double xnum, xden;
  58.  
  59.         if (y > YMAX) {
  60.                 errno = ERANGE;
  61.                 return 0.0;
  62.         }
  63.         if (modf(x*0.63661977236758134308, &xn) >= 0.5)
  64.                 xn += (x < 0.0) ? -1.0 : 1.0;
  65.         f = (x - xn*1.57080078125) + xn*4.454455103380768678308e-6;
  66.         if (fabs(f) < 2.33e-10) {
  67.                 xnum = f;
  68.                 xden = 1.0;
  69.         } else {
  70.                 g = f*f;
  71.                 xnum = P(f,g);
  72.                 xden = Q(g);
  73.         }
  74.         flag |= ((int)xn & 1);
  75.         switch (flag) {
  76.         case 1:         /* A: tan, xn odd */
  77.                 xnum = -xnum;
  78.         case 2:         /* B: cotan, xn even */
  79.                 return xden/xnum;
  80.  
  81.         case 3:         /* C: cotan, xn odd */
  82.                 xnum = -xnum;
  83.         case 0:         /* D: tan, xn even */
  84.                 return xnum/xden;
  85.         }
  86.         return 0.0;
  87. }
  88.  
  89.